############################## bkmr三个分位数的暴露-应答趋势#############################
##### 科研工具5  贝叶斯回归-三个分位数的暴露-应答趋势
library(ggplot2)
library(bkmr)
library(readr)
setwd('父级路径')

demo <- read_csv("source.csv")
demo <- data.frame(demo)
# demoS4 <-demo[,c('Age','Sex','race','edu','FPL','activity','smoke',
#                  'drink','Hypertension','CKD','urinecreatinine','BMI','MECPP',
#                  'MnBP','MEHHP','MEOHP','MiBP','MiNP','MCOP','MCPP','MEP','MBzP',
#                  'Hyperuricemiaint','LBDSUASI')]

tryCatch({
  alldataform
},error = function(e){
  print("获取数据出错")
})

# demoS4 <- demoS4[c(1:200),]
y <- demoS4$Hyperuricemiaint
# Z <- demoS4[,c('MECPP','MnBP','MEHHP','MEOHP','MiBP','MiNP',
#                'MCOP','MCPP','MEP','MBzP')]

tryCatch({
  xdatarep
},error = function(e){
  print("获取数据出错")
})
# X <- demoS4[,c('Age','Sex','race','edu','FPL','activity','smoke',
#                'drink','Hypertension','CKD','urinecreatinine','BMI')]

tryCatch({
  logparseparam
},error = function(e){
  print("log 转换出错提示")
})

set.seed(2012)
fitkm <- kmbayes(y =y, Z = Z, X = NULL, iter = 100, verbose = FALSE, varsel = TRUE,family = 'binomial',est.h = TRUE)
pred.resp.bivar <- PredictorResponseBivar(fit =fitkm , min.plot.dist = 1)
pred.resp.bivar.levels <- PredictorResponseBivarLevels(  
  pred.resp.df = pred.resp.bivar, 
  Z = Z, qs = c(0.1, 0.5, 0.9))


Cairo::CairoTIFF(file="PredictorResponseBivar.tiff", width=8, height=8,units="in",dpi=150)
ggplot(pred.resp.bivar.levels, aes(z1, est)) + 
  geom_smooth(aes(col = quantile), stat = "identity") + 
  facet_grid(variable2 ~ variable1) +
  ggtitle("h(expos1 | quantiles of expos2)") +
  xlab("expos1")
dev.off()